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Abstract 

Numerical issues arising in computations of viscous flows in corners formed by a liquid- 
fluid free surface and a solid boundary are considered. It is shown that on the solid a 
Dirichlet boundary condition, which removes multivaluedness of velocity in the 'moving 
contact-line problem' and gives rise to a logarithmic singularity of pressure, requires a 
certain modification of the standard finite-element method. This modification appears 
to be insufficient above a certain critical value of the corner angle where the numerical 
solution becomes mesh-dependent. As shown, this is due to an eigensolution, which exists 
for all angles and becomes dominant for the supercritical ones. A method of incorporating 
the eigensolution into the numerical method is described that makes numerical results 
mesh-independent again. Some implications of the unavoidable finiteness of the mesh 
size in practical applications of the finite-element method in the context of the present 
problem are discussed. 

1 Introduction 

The ability of a numerical scheme to accurately approximate a physical problem in a domain 
containing corners is critical for the description of a number of phenomena, ranging from 
electromagnetic wave propagation in a waveguide [1] to die-swell effects in polymer extrusion 
p^. Often, an analysis of such problems reveals singular behaviour of variables as the corner 
is approached, which requires special numerical treatment. Such problems are well known and 
have been thoroughly investigated in the setting of fracture mechanics, where one considers the 
propagation of a crack into a material [3], and have also been studied in some fluid dynamics 
problems [4J. 

Our interest here is the viscous flow in a corner formed between a liquid-fluid free surface 
and a solid boundary. Although a free surface is generally bent, to leading order as the corner 
is approached, it is often possible for the purpose of a local analysis to consider the flow domain 
as having a wedge shape. This is the case, for example, in dynamic wetting flows [5], where 
the free surface and the solid boundary form what is referred to as the 'contact angle' and the 
liquid-fluid-solid 'contact line' moves with respect to the solid surface. This differs from the 
situation considered in some previous investigations on flow in a corner where the contact line 
is stationary with respect to the solid and motion is generated by disturbances in the far field 

It is well known that the classical fluid-mechanical approach when applied to dynamic 
wetting problems fails to provide an adequate description of the flow [7j. The conventional 
remedy to the problem is to relax the no-slip boundary condition on the solid surface and allow 
for 'slip' between the liquid and solid. A number of different forms for this slip behaviour have 
been examined in the literature (for a recent review see Ch. 3 of [5J). Broadly, these split into 
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conditions which (i) relate tangential stress to the slip velocity, such as the Navier condition 
[8], or (ii) explicitly prescribe the velocity along the sohd surface. 

In this paper, we consider numerical problems arising in the second case where the velocity 
along the solid surface is a priori prescribed in a form which ensures that a solution exists, and 
that in the far field the usual no-slip condition is restored. This approach is appealing to some 
users due to its mathematical simplicity, and our goal here is to show numerical pitfalls one 
comes across in its numerical implementation and give a method of overcoming them which 
provides a framework for modelling this class of problems. A number of functions prescribing 
the fluid velocity on the solid surface have been proposed in the literature [9l [TOl [11] and here 
we consider just one of these which captures all the main features of the problem. 



2 Problem formulation 

The problem is most easily formulated in a polar coordinate system (r, 6) in a frame moving 
with the contact line (now referred to in our two-dimensional domain as the corner point). The 
wedge is formed by a solid surface at 6' = which moves at speed U parallel to itself, a fiat 
free surface at 6' = a and a 'far field' boundary which is placed at an arc of a sufficiently large 
radius r = R. 

The liquid is Newtonian and incompressible, with density p and viscosity fi. Near the corner 
the flow is characterized by a small length scale so that the Reynolds number Re based on this 
scale is small. Then as Re 0, to leading order in Re we have the Stokes floA^Q The non- 
dimensional Stokes equations for the bulk pressure p and the radial and azimuthal components 
of velocity (m, v) take the form: 

i^ + iS-O. (0<,.<fl.O<.<a.). (1) 
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On the solid surface, for a solution not to have multivalued velocity at the corner point 
[T2| [5], we replace the no-slip condition {u = 1) with a prescribed velocity that has free-slip at 
the corner point and attains no-slip in the far field, that is: 

u = 0, at r = and m — > 1, as r — > cxd. (3) 

Following [11], we use an exponential form for this function and, to complete the boundary 
conditions on the solid surface, it is combined with the usual impermeability condition for the 
component of velocity normal to the surface: 

u=l- exp(-r/s), f = 0, {0 <r < R, 6 = 0). (4) 

The region in which the velocity deviates from no-slip is characterized by the value of s, which 
is a (non-dimensional) 'slip length'. 

On the free surface, we have the standard boundary conditions of zero tangential stress and 
impermeability: 

du 

_=0, ^ = 0, {0<r<R, 9 = a). (5) 

^Thc analysis remains valid for the full Navier-Stokes problem since it considers the r — > limit, and here 
we consider the Stokes equations in the whole region as a convenient way to illustrate the idea. 
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In the far field, we assume that the flow is fully developed and apply 'soft' conditions: 

du dv 

- = - = 0, {r = R,0<e<a), (6) 

which imply that the influence of slip has attenuated; these conditions are satisfied by the 
(multivalued at the corner point) solution obtained using the no-slip condition all along the 
solid surface [13]. 

Equations ([T|-(|6| fully specify the problem of interest. 



3 Local asymptotics 



Consider the leading-order asymptotics for the solution of ([T])-([6]) as r — >• [5] that we will later 
need to use in the numerical code and to provide a test of accuracy of the numerical results 
presented in the next section. After introducing the stream function ijj by 



u 



r 09 ' 



V = — 



dil) 
dr ' 



(7) 



equations ([T])-([2]) are reduced to a biharmonic equation A'^ip = with boundary conditions 
(|4])-([5| taking the form 



r (1 — exp (— r/s)) 



^- = 0, 



0, < r < i?), 



0, ^ = 0, 



a, < r < i?), 



(8) 



(9) 



where, for definiteness, we assign the value zero to the streamline coinciding with the wedge's 
boundary. 

Condition ([S]) is the only inhomogeneous boundary condition in the problem, i.e. the con- 
dition that drives the flow. To leading order as r — 0, it has the form 



dip 



ar 



+ O (r^) 



(10) 



where a = 1/s. An alternative prescribed velocity that satisfies ^ and (10) known in the 
literature [9] is given by u = (r/s)/(l + r/s) and the asymptotic analyses throughout this 



paper are equally valid for this function as well. The form ( 10 ) suggests looking for the leading- 
order term of the local asymptotics in the form ip = r'^F{6), which is a particular case from a 
known family of separable solutions of the biharmonic equation of the form ip = r^F[6). After 
substituting ip = r'^F{6) into A^ip = 0, one arrives at 



^ = (Bi + 320 + ^3 sin 29 + cos 29) , 

where the constants of integration i?j (i = 1, . . . , 4), found from (|8|-(|9]), are given bjj^ 

aa sin 2a 



Bi = — i?4 



B2 = -Bi/a, B-i = Bi cot 2a. 



2a cos 2a — sin 2a 
The pressure field obtained from ^ using ^ and (11) has the form 

p = 4i?2 Inr + Pq. 



(12) 



(13) 



^Here, we correct a typographical error in Bi on p. 126 of [14 and on p. 153 of [5]. 
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where po is a constant which sets the pressure leveL 

It is immediately obvious from (12) that the coefficients are singular when 2a;cos2Q; — 
sin 2a; = 0, which occurs at a critical value Oc determined by tan(2ac) = 2a;c. In the range of 
interest, i.e. for < < 180°, we have ~ 128.7°. 

It is noteworthy, that in the limit r — > 0, the velocity scales linearly with r whilst the 
pressure is logarithmically singular at the corner and is independent of the angular coordinate 
9. 



4 Numerical results 

From a numerical viewpoint, the steady fixed-boundary problem considered in this paper is 



complicated only by the presence of a singularity in the pressure which, according to (13), is 
logarithmic as r — 0. The simplicity of the rest of the problem and the availability of asymptotic 
results, which not only give the behaviour of the velocity and pressure near the corner, but also 
provide the coefficients, make this a perfect testing ground for a numerical method's ability to 
approximate flows in corner regions formed by boundaries on which different types of boundary 
conditions are applied. 

In the standard implementation of the finite-element, as well as finite-difference, algorithm, 
one assigns an a priori unknown finite value to the pressure at the corner point. If such a code 
attempts to approximate a solution where the pressure at the corner point is singular, like the 
one whose local asymptotics we considered earlier, the nodal value of the pressure, as well as 
the pressure at the neighbouring nodes, will vary as one refines the mesh. In other words, an 
attempt to approximate a singular analytic solution using regular numerical representations of 
the unknown function on each element will lead to a numerical 'solution' that is mesh-dependent 
and hence, strictly speaking, it is not a solution to the original problem formulated in terms of 
PDEs that 'do not know' about any mesh. 

In order to achieve a uniformly valid solution in the framework of a finite-element method, 
one approach is to use singular elements, i.e. to redefine the pressure interpolation in the 
elements that contain the corner point node in such a way that the pressure is allowed to 
be infinite at the corner point and behave as described by the local asymptotics. A simple 
implementation of this idea is described in 

In the present work, the problem formulated in Section [2] has been considered using part 
of a finite-element-based numerical platform which has been developed to simulate a range of 
microfluidic capillary flows and has already been used to obtain new results for the flow of 
liquids over surfaces of varying wettability fTEl ITU] . The idea here, is to modify the finite- 
element's basis function $p associated with the pressure at the corner point. In the standard 
triangular Taylor-Hood element with six velocity nodes and three pressure nodes, $p is linear 
and takes the value 1 at the corner point and at all other pressure nodes. Now, instead of 
using $p and determining the coefficient in front of it, we will be using and determining the 
coefficient in front of a singular basis function $5 = $p Inr. We denote the computed coefficient 
of $s as Bp. Instead of $p one could use other functions to pre-multiply the logarithm, for 
example sin('7r$p/2); such functions have been tried and it was found that they do not offer 
noticeable advantages over $p. 

It is pointed out in |1] that the usual Gaussian quadrature is not well suited to the inte- 
gration of a singular function, such as Inr, and a special Simpson quadrature routine has been 
suggested to provide a very accurate approximation of the integrals. This routine has also been 
incorporated into our numerical platform; however, we found that using Gaussian quadrature 
with enough integration points provided an accurate enough estimation of the integral and was 
significantly quicker and easier to implement. Sixteen Gauss integration points were found to 
be more than sufficient. 



4 



0=75 




Figure 1: Left: streamlines in increments oi ifj = 0.05, with the 6 = 0,a interfaces correspond- 
ing to = 0. Right: comparison of the computed velocity with the analytic prediction (dashed 
line) along the liquid-solid (1) and liquid-fluid (2) interfaces. 

At the critical angle ac ~ 127.8°, more complex asymptotic analysis should be considered 
to resolve the singular behaviour and the results should be incorporated into the code in a 
way similar to what is being described. However, since this paper is concerned with the gen- 
eral numerical approximation of corner flows which contain singularities and their numerical 
treatment, we shall consider angles away from this critical value, i.e. the range where, as one 
might expect at this stage, our asymptotics of Section |3] can be used. First, we shall consider 
subcritical angles a < ac, in particular a = 75° as a representative case. We take s = 0.1 and 
i? = 10 in all the simulations that we present as an investigation into the variation of these 
parameters would be about the physical problem rather than its numerics and would detract 
from the main emphasis of this paper. 

4.1 Numerical approximation at subcritical corner angles 

In Fig. [T| we show the streamlines generated by the exponential slip model. As expected, the 
prescribed velocity on the solid draws fluid near the solid out of the corner, thus reducing the 
pressure there which then sucks in fluid from the far field. In the same figure, the components 
of the radial velocity along the interfaces are compared to the asymptotic predictions; the 
agreement between the calculated velocity along the liquid-fluid interface and the asymptotic 
prediction is visibly excellent (agreement along the liquid-solid interface near the corner is 
guaranteed as the velocity is prescribed, so we give it here just to show the range in which the 
velocity varies). 

The plot of pressure along the two interfaces in Fig. [2] shows that the pressure is indeed 
6'-independent, and it is almost graphically indistinguishable from the asymptotic result. To 
confirm the mesh- independence of the result, we also consider how well the singular behaviour 
of pressure is captured as the mesh is resolved over ten orders of magnitude using ten different 
meshes, characterized by the width of the smallest element ri. To do so, we show both the 
coefficient Bp in front of the basis function $s (curve 1 in Fig. |2]) and the appropriate gradient 
determined from the pressures pi, p2 at the two pressure nodes on the solid surface, for simplic- 
ity, closest to the corner point at radial distances ri, r2, which is given by {p2—pi)/{lnr2 — \nri) 
(curve Ig). This second method allows us to compare the code with the singular elements to 
one without (curve 2g in the plot) where Bp is not explicitly calculated. 

The clearest conclusions are drawn from the results of the second method (curves Ig and 
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Figure 2: Left: comparison of computed pressure along the 6 = 0,a interfaces, which are 
graphically indistinguishable, as the corner is approached with the analytic pressure (dashed 
line). Right: convergence of the pressure coefficient of Inr with (1) and without (2) singular 
elements as the mesh is resolved. Bp is calculated from the singular element's coefficient (1) 
and from the local gradient (Ig) and (2g). The analytic value is 7.23 and ri is the size of the 
smallest element. 
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Figure 3: Left: streamlines in increments of = 0.1 with the 6 = 0,a interfaces corresponding 
to ip = 0. Right: comparison of the velocity along the liquid-fluid interface with the analytic 
prediction (dashed line). 

2g). Here, the plot shows that, by using the singular element, the local gradient converges to the 
asymptotic value of Bp = 7.23 (the dashed line in the figure): without these singular elements 
the code converges to the incorrect value as the mesh is resolved. This is strong support for 
the inclusion of singular elements in dynamic wetting codes. Without them, the behaviour of 
pressure is wrong not only in the element adjacent to the corner point, but, by continuity, also 
in a neighbourhood of this element. 

For a < ac the asymptotics and numerics are in excellent agreement, and the special 
treatment of the corner singularity was a success. Now we consider a supercritical angle a = 
175° > ac- 

4.2 Numerical approximation at supercritical corner angles 

The streamlines in Fig. |3] are as one may intuitively expect, but when we compare the numerical 
and asymptotic results for the velocity along the liquid-fluid interface there is no agreement. In 
fact, the asymptotic result predicts that the flow should be up the liquid-fluid interface, which 
is clearly not the case in the computed solution. 
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Figure 4: Left: computed pressure in the vicinity of the corner along the liquid- solid (1) and 
liquid-fluid (2) interfaces compared to the analytic prediction (dashed line). Right: pressure 
coefficient compared to the analytic value 1.122 (dashed line) plotted against ri the size of the 
smallest element. 



The computed pressure along the liquid-solid interface is given as curve 1 in Fig. |4j It is 
not only that the pressure strongly deviates from the analytic result (dashed line) as the corner 
point is approached; one can see that there also appear huge oscillations: the pressure decrease 
in the element adjacent to the special corner elements (the line to the left of the point 10~^ in 
the plot) is followed by a steep increase in the element comprising the corner point (not shown 
in the semilogarithmic plot). Such mesh-dependence of the numerical result indicates that the 
obtained solution cannot be regarded as a valid approximation of the solution to the original 
set of partial differential equations. This conclusion is re-enforced when we study the value of 
Bp, the coefficient to the logarithm in the singular element, as we refine the mesh: there is no 
convergence. A similar trend is observed if we study Bp using the local gradient method. The 
same conclusions may be drawn for all supercritical angles. 

Thus, the standard FEM coupled with the local-asymptotics-based approximation of the 
pressure does not allow one to obtain an acceptable numerical approximation for solutions 
of the Stokes equations in a corner with a combination of Dirichlet and Neumann boundary 
conditions on the interfaces for all angles greater than 127.8°. 

The robustness of the obtained numerical 'solution' suggests that there is a fundamental 
numerical problem. It should be emphasized that we arrived at this difficulty just by varying 
the wedge angle in a code that produces excellent results for smaller wedge angles, and then 
cannot provide any mesh-independent solution after a critical angle. An immediate (tentative) 
explanation for this situation is that, besides the analytical solution whose asymptotics has 
been considered in Section [3] and incorporated into the code, there exists a 'local eigensolution', 
i.e. a solution satisfying zero boundary conditions on the sides of the wedge, that becomes 
dominant for the supercritical angles. We will now examine this conjecture. 

4.3 Asymptotics of an eigensolution 

The near- field asymptotics of the eigensolution to the biharmonic equation satisfying conditions 

dip d'^ip 
= ^ = 0, on 9 = 0, and ip = —r^r = 0, on 9 = a, (14) 
o9 o9^ 

is given by 

ipe = [Ai sin (A^) + A2 cos {X9) + A3 sin ((A - 2) ^) + A4 cos ((A - 2) 9)] . (15) 
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Using (14) and noting that A 7^ 1, 2, we find that A is determined by the equation: 



2 sin (Aa) cos ((A - 2) a) = A sin (2a) . (16) 



Defining the degree of freedom hj A = Ai, the boundary conditions (14) give: 



^ A(A-2)-^sin((A-2)a)-sin(Aa) A 
- "^^ - cos ((A -2) a)- cos (Aa) ' " "^A^' ^^^^ 

and the pressure has the form: 

Pe = 4A(A - ly-^ [as cos((A - 2)6) - sin((A - 2)^)] + pi (18) 
where = Ai/A, i = 3,4 and pi is a constant setting the pressure leveL Promisingly, equation 



([16]) has roots A G (1, 2) for a e (a^, 180°), with A ^ 2 as a ^ a^, A ^ 3/2 as a 180° and 
A as a function of a varying monotonically between these hmiting values. 

This eigensolution has been derived in a number of other works considering the flow in 
a corner formed between an impermeable no-slip boundary and an impermeable, sometimes 
free, zero-tangential stress boundary, e.g. in [17]; these are sometimes referred to as 'stick-slip 
phenomena' flEl |6]. The difference between our problem and the aforementioned flows with 
static corner points is that, unlike these situations where the flow is driven by the far field, 
here the fluid motion is generated by the movement of the solid and analytically this behaviour 
is captured in the asymptotics of Section |3| The eigensolution comes on top of this solution 
and, in the near fleld, it 'does not know' about the motion of the solid, although, ultimately, 
it is the solid's motion that generates the flow in the far fleld that gives rise to this solution. 
The eigensolution exists in the range of subcritical angles as well, but there it is regular in all 
variables and therefore causes no problem for numerical computations; it is only for a > Oc 
that the eigensolution becomes both singular and dominant. 

For a < ttc the pressure at the corner point can be referred to as single- valued: the coefficient 
in front of the logarithm is independent of 6. In contrast, the solution for a > etc is manifestly 



multivalued as predicted by the eigensolution (18) and as seen numerically: if one takes a 
vicinity of the corner point, then, no matter how small this vicinity is, there will be points 
which are equidistant from the corner point with an arbitrarily large pressure difference. Here, 
being interested in the numerical side of the problem, we set aside physical arguments that 
might arise in connection with the obtained solution. 

The existence of an eigensolution and its dominance for a > ac suggests that our initial 
attempt at computing the flow at large angles were flawed because, given the asymptotics 
of Section [3} we assumed that the pressure scaled as Inr whereas in fact the most singular 
term (i) has order r"*' where A; = 2 — A G (0, 1) and (ii) is dependent on 6. This suggests a 
generalisation of our approach: we need to incorporate the new singular behaviour into the 
special elements adjacent to the corner point. Due to the presence of the eigensolution, we 
now have an unknown constant A in our asymptotics which will prevent us from comparing a 
priori determined analytic curves with our numerical results. However, once A is determined 
numerically, we may use it to extrapolate the asymptotic behaviour outside the singular element 
in which it is calculated, i.e. we may then compare, a now semi-analytic, asymptotic prediction 
to the computed solution globally. 

An alternative method that we used to verify our singular element solution and do not 
describe in detail here is to analytically remove the eigensolution, which is the cause of numer- 
ical difficulties, prior to computation and then superimpose it back on after. This approach 
that has been shown to be essential for the simulation of flows using the Navier slip condition, 
a Robin- type boundary condition, on the solid surface is described elsewhere [19]. In more 
complex problems, where the corner is just one element, this method of removing the eigen- 
solution everywhere is overly complex and a local method should be used which removes the 
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Figure 5: Left: computed pressure distributions in the vicinity of the corner along the hquid- 
sohd (1) and hquid-fluid (2) interfaces compared with the semi-analytic result (dashed lines). 
Right: convergence of pressure coefficient Ap to 0.092, plotted against ri the size of the smallest 
element. 

eigensolution near to the corner; this method is also described in ^19j using examples of full 
scale dynamic wetting simulations. 



4.3.1 Modified singular elements 



In the limit as r — 0, the first two terms in the expansion of pressure are p = Apr''^~'^g{6) + 
Bplnr. For a > Oc both terms are singular. We will begin by using only the leading order 
term in r in our singular elements and will return to the two-term approximation later. Taking 
the leading term, our new singular elements have a basis function of the form: 



$p In r. 



a < Or 



A-2 



g{6), a> ac 



(19) 



where the unknown coefficients are Bp and Ap, respectively. 
In Fig. Is] with a = 175° and hence, from (16), A 



1.529, we see that the implementation 
of the new singular elements solves our previous problems by (i) removing the oscillations in 
pressure as the corner point is approached, and (ii) converging as the mesh is refined. 

The value of A = Ap, which is the coefficient of the singular basis functions and is deter- 
mined by the finite element method, may be used after computation with the supplementary 



asymptotics of Section |4.3| to produce a fully determined asymptotic solution for the velocity 
and pressure. Then we may compare the analytic results of Section 4.3 with our numerical 



results globally. It should be pointed out that using the value of A to extrapolate the analytic 
behaviour of the eigensolution well outside the first elements provides a quick check to see if 
A is in the correct range, this value does not in any way actually determine the velocity field 
outside the first elements. The comparison with pressure in Fig. [5] shows good agreement be- 
tween numerics and asymptotics; however, very close to the corner point along the free surface 
the numerical solution is not as smooth as the asymptotic result. This is no surprise given the 
huge gradients in pressure which are being approximated by linear basis functions both in the 
radial and angular directions. 

In Fig. [6} we compare the velocity along the interfaces of the computed numerical solution 
to the asymptotic result, showing in particular how the full asymptotic solution is a super- 
position of the eigensolution Ue and the supplementary solution Ua- We see that, although 
the supplementary asymptotics predicts that fiow will be reversed near the contact line, this 
is blown away by the strength of the eigensolution, which restores what one would intuitively 
think is the correct direction for the fiow. The agreement we see in this figure is sufficient when 
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Figure 6: Comparison of the computed velocity along the liquid-fluid interface u with the 
semi-analytic result which is shown decomposed into a contribution from the supplementary 
asymptotics Ua and the eigensolution contribution Ue with A=0.092. 

we consider that it is determined by the coefficient of the singular pressure which only plays a 
role in elements adjacent to the corner point. 

4.3.2 Numerics incorporating two-term asymptotics 

For an angle of a = 175° we have completely resolved the situation. However, we have observed 
for smaller angles, roughly 150°, that, in terms of mesh refinement, convergence is 

very slow, i.e. a mesh-independent regime is only realised for exceptionally well resolved meshes 
which are well outside the scope of most numerical platforms where the corner is just one part 
of a larger problem. In this range of angles, although asymptotically in the limit r — »• 
the logarithmic pressure behaviour (p ~ Inr) is overshadowed by the eigensolution, where 
p r^~'^g{6), in reality, unavoidable finiteness of the resolution of the mesh means that one 
could be sufficiently far away from the corner point for the logarithm to be still dominant in 
the numerics. To see if this is indeed the case, we consider the relative size of the two pressure 
terms at the edge opposite the corner point in the first elements. Taking this element to have 
size 10~", we see that the two singular functions are comparable, at a critical value Ac, at the 
edge of the first element when 

InlO"" = 10-"(^=-2), (20) 

that is when 

In I — nln lOl , , 

Ac = 2 -. (21) 

nlnlO ^ ^ 



Taking n = 6, from (21) we have that Ac = 1.81, which means, using (16), that the logarithmic 
behaviour dominates a numerical scheme, with the stated spatial resolution, in the range of 
angles etc < « < 143°. Thus, we have seen both numerically and analytically, that even for an 
extremely well resolved mesh, we should expect that there is a range of angles in which the 
logarithmic solution cannot be neglected; approximately for a < 150°. This is a very serious 
issue which must be resolved as almost all numerical schemes will not be able to afford the 
required resolution to accurately capture the singular behaviour. 

The solution to this issue is to include the second term in the asymptotic expansion of 
pressure, so that p = Ap^pr''^~'^g{6) + Bp^plnr for a > ac- The coefficient to this logarithm 
Bp is known exactly from the supplementary asymptotics (13) and we have already shown in 



Section |4.1| that this value is numerically reproduced. Therefore, the simplest solution is to 
prescribe the value of Bp and see if this improves the speed of convergence of Ap. 

In Fig.[7]we show how nicely this method works with, as one would expect, the most improve- 
ment occurring for smaller angles where the logarithmic pressure is strongest. For example, for 
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Figure 7: A comparison of the convergence of the pressure coefficient Ap with mesh size ri. 
With (dashed hue) and without (full line) the supplementary asymptotics prescribed, for a 
range of different angles. 

a = 150°, where the converged value is Ap = 1.007, if we neglect the logarithm, then we require 
~ 5000 elements to get within 5% of this value, but, if we include the asymptotic logarithmic 
behaviour then we only need ~ 2500 elements to attain the same accuracy. Computationally 
this is a terrific saving and, given the simplicity of its implementation, we conclude that the 
logarithmic asymptotic behaviour should be included for all angles a > ac- 

5 Conclusion 

We have shown that accurate numerical calculation of viscous flows near corners of the flow 
domain in the framework of the finite-element method requires special treatment of the elements 
adjacent to the corners. In the case of zero-stress/prescribed velocity boundary conditions on 
the sides of the corner, the range of corner angles is spUt into two distinct regions. For the 
angles below the critical angle ~ 128.7, it is sufficient to use a logarithmic basis function 
for the pressure, whereas for supercritical angles there appears a 'hidden' eigensolution which 
considerably complicates numerics. In order to obtain an acceptable (i.e. mesh-independent) 
solution, one has to incorporate this eigensolution into the code by altering the basis functions 
for the pressure in the elements adjacent to the corner point. Close to the critical angle it 
becomes necessary to use a two-term asymptotics in the numerical algorithm, including the 
leading terms of both the eigensolution and the solution of the inhomogeneous problem, as the 
finiteness of the mesh size could result in the algebraically and logarithmically singular terms 
having comparable values. 

An issue that is now opened up for numerical and analytic investigation is how to generalize 
the developed methods for a three-dimensional case, i.e. in a situation where both the contact 
angle and the direction of velocity of the solid vary along a contact line. The first of these aspects 
becomes particularly challenging when the angle varies from subcritical to supercritical. 
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